/*================================================================= *
 * GDALINFO.C
 *     Mexfile that returns metadata about the gdal file.
 *
 * USAGE:
 *
 *    METADATA = MEXGDAL(GDALFILE);
 *
 * Output:
 *
 * In case of an error, an exception is thrown.
 *
 *=================================================================*/
/* $Revision: 1.4 $ */
#include "gdal.h"

#include "mex.h"

int record_geotransform ( char *gdal_filename, GDALDatasetH hDataset, double *adfGeoTransform );
void handle_overviews ( GDALRasterBandH hBand, mxArray *band_struct );
int unpack_verbose ( const mxArray *field );
mxArray * populate_metadata_struct ( char * );

/*
 * If this flag is tripped, then we want to provide debugging output.
 */
int mexgdal_verbose;

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )
{ 
    /*
     * indices
     * */
    int i, j, index, t_index;

    /*
     * Construct error and warning messages using this buffer.
     * */
    char error_msg[500];


    /*
     * name of input file
     * */
    char *gdal_filename;

    /*
     * short cut to the mxArray data
     * */
    double *dptr;

    /*
     * Pointer to matlab raster arrayy
     * */
    mxArray *mxGDALraster;

    /*
     * Pointers to matlab array aliases.
     * */
    mxArray *mx_input_gdal_file;


    int status;            /* success or failure */

    /*
     * pointer structure used to query the gdal file.
     * */
    GDALDatasetH hDataset;


    GDALRasterBandH hBand;

    /*
     * GDT Byte?, GDT UInt32?  What is it?
     */
    GDALDataType gdal_type;

    /*
     * Which band do we want to retrieve?
     * */
    int requested_band;

    /*
     * What overview are we to retrieve?  If any at all?
     * */
    int requested_overview;

    /*
     * How big is the output datatype?  We can have either Byte or Float64 for now.
     * This is in bytes.
     */
    int out_type_size = 0;
    unsigned char  *transpose_buffer_byte;
    double  *transpose_buffer_float64;

    /*
     * The data is to be copied into void_buffer initially.  We then transpose it
     * otherwise it comes out, well, transposed.  
     */
    void   *void_buffer;
    void   *transpose_void_buffer;

    /*
     * size of allocated matlab array.
     */
    int rasterDims[2];

    /*
     * This flag keeps track of whether the default assumptions about the
     * output should be followed.  If no 2nd input argument is given, then
     * this is the case.
     */
    int defaults_are_invoked;

    /*
     * If this flag is tripped, then we only want to return the metadata.
     * */
    int gdal_dump;

	/*
	 * Length of string.
	 */
	size_t buflen;


    /*
     * Set up the defaults.
     */
    defaults_are_invoked = 0;  /* Assume the user is going to provide input options. */
    gdal_dump = 0;             /* We aren't looking for metadata only. */
    requested_band = 1;        /* Get the first band unless we are told otherwise. */
    requested_overview = -1;   /* Don't get any overview unless specifically asked for. */
    mexgdal_verbose = 1;       /* Don't provide debugging output unless told otherwise. */

    /* 
     * Check for proper number of arguments
     */
    if (nlhs != 1) { 
        mexErrMsgTxt("Only one output argument, please."); 
    }
    if (nrhs != 1) { 
        mexErrMsgTxt("Only one input argument, please."); 
    }

    /*
     * The first argument must be character.
     * */
    if ( mxIsChar(prhs[0]) != 1 ) {
        mexErrMsgTxt("Input file name must be a string\n" );
    }

    buflen = mxGetNumberOfElements(prhs[0]) + 1; 
    gdal_filename = (char *) mxCalloc ( buflen, sizeof(char) );



    /*
     * copy the string data from prhs[0] into a C string.
     * */
    status = mxGetString(prhs[0],gdal_filename,buflen);
    if ( status != 0 ) {
        mexErrMsgTxt ( "Not enough space for input file argument.\n" );
    }


    GDALAllRegister();

    plhs[0] = populate_metadata_struct ( gdal_filename );
    return;



}



/*
 * POPULATE_METADATA_STRUCT
 *
 * This routine just queries the GDAL raster file for all the metadata
 * that can be squeezed out of it.
 *
 * The resulting matlab structure is by necessity nested.  Each raster 
 * file can have several bands, e.g. PNG files usually have 3, a red, a 
 * blue, and a green channel.  Each band can have several overviews (tiffs
 * come to mind here).
 *
 * Fields:
 *    ProjectionRef:  a string describing the projection.  Not parsed.
 *    GeoTransform:  
 *        a 6-tuple.  Entries are as follows.
 *            [0] --> top left x
 *            [1] --> w-e pixel resolution
 *            [2] --> rotation, 0 if image is "north up"
 *            [3] --> top left y
 *            [4] --> rotation, 0 if image is "north up"
 *            [5] --> n-s pixel resolution
 *
 *    DriverShortName:  describes the driver used to query *this* raster file
 *    DriverLongName:  describes the driver used to query *this* raster file
 *    RasterXSize, RasterYSize:  
 *        These are the primary dimensions of the raster.  See "Overview", though.
 *    RasterCount:
 *        Number of raster bands present in the file.
 *    Driver:
 *        This itself is a structure array.  Each element describes a driver
 *        that the locally compiled GDAL library has available.  So you recompile
 *        GDAL with new format support, this structure will change.
 *
 *        Fields:
 *            DriverShortName, DriverLongName:  
 *                Same as fields in top level structure with same name.
 *
 *    Band:
 *        Also a structure array.  One element for each raster band present in 
 *        the GDAL file.  See "RasterCount".
 *
 *        Fields:
 *            XSize, YSize:
 *                Dimensions of the current raster band.
 *            Overview:
 *                A structure array, one element for each overview present.  If
 *                empty, then there are no overviews.
 *            NoDataValue:
 *                When passed back to MATLAB, one can set pixels with this value
 *                to NaN.
 *
 * */
mxArray * populate_metadata_struct ( char *gdal_filename ) {
    /*
     * Number of available drivers for the version of GDAL we are using.
     * */
    int driverCount;

    /*
     * These are used to define the metadata structure about available GDAL drivers.
     * */
    char *driver_fieldnames[100];  
    int num_driver_fields;

    mxArray *driver_struct;
    mxArray *mxtmp;
    mxArray *mxProjectionRef;
    mxArray *mxGeoTransform;
    mxArray *mxGDALDriverShortName;
    mxArray *mxGDALDriverLongName;
    mxArray *mxGDALRasterCount;
    mxArray *mxGDALRasterXSize;
    mxArray *mxGDALRasterYSize;

    /*
     * These will be matlab structures that hold the metadata.
     * "metadata_struct" actually encompasses "band_struct", 
     * which encompasses "overview_struct"
     * */
    mxArray *metadata_struct;
    mxArray *band_struct;

    /*
     * Loop indices
     * */
    int j, band_number;

    /*
     * short cut to the mxArray data
     * */
    double *dptr;

    /*
     * This is the driver chosen by the GDAL library
     * to query the dataset.  
     * */
    GDALDriverH hDriver;

    /*
     * pointer structure used to query the gdal file.
     * */
    GDALDatasetH hDataset;
    GDALRasterBandH hBand;

    /*
     * Construct error and warning messages using this buffer.
     * */
    char error_msg[500];


    int status;            /* success or failure */

    /*
     * bounds on the dataset
     * */
    double     adfGeoTransform[6];

    /*
     * number of fields in the metadata structures.
     * */
    int num_struct_fields;   
    int num_band_fields;

    /*
     * this array contains the names of the fields of the metadata structure.
     * */
    char *fieldnames[100];  
    char *band_fieldnames[100];  

    /*
     * Dimensions of the dataset
     * */
    int xSize, ySize, raster_count;

    /*
     * Datatype of the bands.
     * */
    GDALDataType gdal_type;

    /* temporary value */
    double tmpdble; 

    /*
     * Retrieve information on all drivers.
     * */
    driverCount = GDALGetDriverCount();

    /*
     * Open the file.
     * */
    hDataset = GDALOpen ( gdal_filename, GA_ReadOnly );
    if ( hDataset == NULL ) {
        sprintf ( error_msg, "Unable to open %s.\n", gdal_filename );
        mexErrMsgTxt ( error_msg );
    }

    /*
     * Create the metadata structure 
     * Just one element, with XXX fields.
     * */
    num_struct_fields = 9;
    fieldnames[0] = strdup ( "ProjectionRef" );
    fieldnames[1] = strdup ( "GeoTransform" );
    fieldnames[2] = strdup ( "DriverShortName" );
    fieldnames[3] = strdup ( "DriverLongName" );
    fieldnames[4] = strdup ( "RasterXSize" );
    fieldnames[5] = strdup ( "RasterYSize" );
    fieldnames[6] = strdup ( "RasterCount" );
    fieldnames[7] = strdup ( "Driver" );
    fieldnames[8] = strdup ( "Band" );
    metadata_struct = mxCreateStructMatrix ( 1, 1, num_struct_fields, (const char **)fieldnames );

    num_driver_fields = 2;
    driver_fieldnames[0] = strdup ( "DriverLongName" );
    driver_fieldnames[1] = strdup ( "DriverShortName" );
    driver_struct = mxCreateStructMatrix ( driverCount, 1, num_driver_fields, (const char **)driver_fieldnames );
    for ( j = 0; j < driverCount; ++j ) {
        hDriver = GDALGetDriver( j );
        mxtmp = mxCreateString ( GDALGetDriverLongName(hDriver) );
        mxSetField ( driver_struct, j, (const char *) "DriverLongName", mxtmp );

        mxtmp = mxCreateString ( GDALGetDriverShortName(hDriver) );
        mxSetField ( driver_struct, j, (const char *) "DriverShortName", mxtmp );

    }
    mxSetField ( metadata_struct, 0, "Driver", driver_struct );

    /*
     * Record the ProjectionRef.
     * */
    mxProjectionRef = mxCreateString ( GDALGetProjectionRef( hDataset ) );
    mxSetField ( metadata_struct, 0, "ProjectionRef", mxProjectionRef );




    /*
     * Record the geotransform.
     * */
    mxGeoTransform = mxCreateNumericMatrix ( 6, 1, mxDOUBLE_CLASS, mxREAL );
    dptr = mxGetPr ( mxGeoTransform );

    status = record_geotransform ( gdal_filename, hDataset, adfGeoTransform );
    if ( status == 0 ) {
        dptr[0] = adfGeoTransform[0];
        dptr[1] = adfGeoTransform[1];
        dptr[2] = adfGeoTransform[2];
        dptr[3] = adfGeoTransform[3];
        dptr[4] = adfGeoTransform[4];
        dptr[5] = adfGeoTransform[5];
        mxSetField ( metadata_struct, 0, "GeoTransform", mxGeoTransform );
    } else {
        sprintf ( error_msg, "No internal georeferencing exists for %s, and could not find a suitable world file either.\n", gdal_filename );    
        mexWarnMsgTxt ( error_msg );
    }

    /*
     * Get driver information
     * */
    hDriver = GDALGetDatasetDriver( hDataset );

    /*
    printf( "Driver: %s/%s\n",
        GDALGetDriverShortName( hDriver ),
        GDALGetDriverLongName( hDriver ) );
        */
    mxGDALDriverShortName = mxCreateString ( GDALGetDriverShortName( hDriver ) );
    mxSetField ( metadata_struct, 0, (const char *) "DriverShortName", mxGDALDriverShortName );

    mxGDALDriverLongName = mxCreateString ( GDALGetDriverLongName( hDriver ) );
    mxSetField ( metadata_struct, 0, (const char *) "DriverLongName", mxGDALDriverLongName );

    xSize = GDALGetRasterXSize( hDataset );
    mxGDALRasterXSize = mxCreateDoubleScalar ( (double) xSize );
    mxSetField ( metadata_struct, 0, (const char *) "RasterXSize", mxGDALRasterXSize );

    ySize = GDALGetRasterYSize( hDataset );
    mxGDALRasterYSize = mxCreateDoubleScalar ( (double) ySize );
    mxSetField ( metadata_struct, 0, (const char *) "RasterYSize", mxGDALRasterYSize );

    raster_count = GDALGetRasterCount( hDataset );
    mxGDALRasterCount = mxCreateDoubleScalar ( (double)raster_count );
    mxSetField ( metadata_struct, 0, (const char *) "RasterCount", mxGDALRasterCount );
    /*
    printf( "Size is %dx%dx%d\n", xSize, ySize, raster_count );
    */

    /*
     * Get the metadata for each band.
     * */
    num_band_fields = 5;
    band_fieldnames[0] = strdup ( "XSize" );
    band_fieldnames[1] = strdup ( "YSize" );
    band_fieldnames[2] = strdup ( "Overview" );
    band_fieldnames[3] = strdup ( "NoDataValue" );
    band_fieldnames[4] = strdup ( "DataType" );
    band_struct = mxCreateStructMatrix ( raster_count, 1, num_band_fields, (const char **)band_fieldnames );

    for ( band_number = 1; band_number <= raster_count; ++band_number ) {

        hBand = GDALGetRasterBand( hDataset, band_number );

        mxtmp = mxCreateDoubleScalar ( (double) GDALGetRasterBandXSize( hBand ) );
        mxSetField ( band_struct, 0, "XSize", mxtmp );

        mxtmp = mxCreateDoubleScalar ( (double) GDALGetRasterBandYSize( hBand ) );
        mxSetField ( band_struct, 0, "YSize", mxtmp );
    
        gdal_type = GDALGetRasterDataType ( hBand );

        mxtmp = mxCreateString ( GDALGetDataTypeName ( gdal_type ) );
        mxSetField ( band_struct, 0, (const char *) "DataType", mxtmp );

        tmpdble = GDALGetRasterNoDataValue ( hBand, &status );
        mxtmp = mxCreateDoubleScalar ( (double) (GDALGetRasterNoDataValue ( hBand, &status ) ) );
        mxSetField ( band_struct, 0, "NoDataValue", mxtmp );
    
        /*
         * Can have multiple overviews per band.
         * */
        handle_overviews ( hBand, band_struct );


    }

    mxSetField ( metadata_struct, 0, "Band", band_struct );

    GDALClose ( hDataset );

    return ( metadata_struct );
}










/*
 * HANDLE_OVERVIEWS
 *
 * If the raster file has overviews, then we need to populate the 
 * metadata structure appropriately.
 * */
void handle_overviews ( GDALRasterBandH hBand, mxArray *band_struct ) {

    /*
     * Number of metadata items for each overview structure.
     * */
    int num_overview_fields = 2;

    char *overview_fieldnames[2];  

    /*
     * Just a temporary matlab array structure.  We don't keep it around.
     * It will hold the X and Y sizes of the overviews.
     * */
    mxArray *mxtmp;

    /*
     * This matlab structure will contain the size of the overview
     * in fields "XSize" and "YSize"
     * */
    mxArray *overview_struct;



    /*
     * GDAL datastructure for the overview.
     * */
    GDALRasterBandH overview_hBand;

    /*
     * Dimensions of the overview
     * */
    int xSize, ySize;


    /*
     * Number of overviews in the current band.
     * */
    int num_overviews;


    /*
     * loop index
     * */
    int overview;


    /*
     * These are the only fields defined for the overview metadata.
     * */
    overview_fieldnames[0] = strdup ( "XSize" );
    overview_fieldnames[1] = strdup ( "YSize" );


    num_overviews = GDALGetOverviewCount( hBand );
    if ( num_overviews > 0 ) {

        overview_struct = mxCreateStructMatrix ( num_overviews, 1, num_overview_fields, 
                    (const char **)overview_fieldnames );

        for ( overview = 0; overview < num_overviews; ++overview ) {
            overview_hBand = GDALGetOverview ( hBand, overview );

            xSize = GDALGetRasterBandXSize ( overview_hBand );
            mxtmp = mxCreateDoubleScalar ( xSize );
            mxSetField ( overview_struct, overview, "XSize", mxtmp );

            ySize = GDALGetRasterBandYSize ( overview_hBand );
            mxtmp = mxCreateDoubleScalar ( ySize );
            mxSetField ( overview_struct, overview, "YSize", mxtmp );
        }
        mxSetField ( band_struct, 0, "Overview", overview_struct );
    }

}




/*
 * record_geotransform:  
 *
 * If the gdal file is not internally georeferenced, try to get the world file.
 * Returns -1 in case no world file is found.
 * */
int record_geotransform ( char *gdal_filename, GDALDatasetH hDataset, double *adfGeoTransform ) {

	int status = -1;
	char generic_buffer[5000];

	if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None ) {
		/*
		printf( "Origin = (%.6f,%.6f)\n",
			adfGeoTransform[0], adfGeoTransform[3] );
		printf( "Pixel Size = (%.6f,%.6f)\n", adfGeoTransform[1], adfGeoTransform[5] );
		*/
		return ( 0 );
	}

	/*
	 * Try a world file.  First the generic extension.
	 * If the gdal_filename is, say, "a.tif", then this
	 * will look for "a.wld".  
	 * */
	if ( GDALReadWorldFile ( gdal_filename, "wld", adfGeoTransform ) ) {
		return ( 0 );
	} 

	/*
	 * Try again, but try "a.tif.wld" instead.
	 * */
	sprintf ( generic_buffer, "%s.xxx", gdal_filename );
	status = GDALReadWorldFile ( generic_buffer, "wld", adfGeoTransform );
	if ( status == 1 ) {
		return ( 0 );
	}

	/*
	 * Newer versions of GDAL will try to guess if you pass NULL.  Older
	 * versions with barf, so becareful about attempting this.
	 * */
	if ( GDAL_VERSION_NUM >= 1210 ) {
		if ( GDALReadWorldFile ( gdal_filename, NULL, adfGeoTransform ) ) {
			return ( 0 );
		}

	}
	return ( -1 );
}

